A monotone nonlinear finite volume method for advection–diffusion equations on unstructured polyhedral meshes in 3D

We present a new monotone finite volume method for the advection–diffusion equation with a full anisotropic discontinuous diffusion tensor and a discontinuous advection field on 3D conformal polyhedral meshes. The proposed method is based on a nonlinear flux approximation both for diffusive and advective fluxes and guarantees solution non-negativity. The approximation of the diffusive flux uses the nonlinear two-point stencil described in [9]. Approximation of the advective flux is based on the second-order upwind method with a specially designed minimal nonlinear correction [26]. The second-order convergence rate and monotonicity are verified with numerical experiments. The discrete maximum principle (DMP) and local mass conservation are important properties of a numerical scheme for the approximate solution of the steady state advection–diffusion equation. An accurate discretization method satisfying DMP is hard to develop. We address the monotonicity condition as the simplified version of the DMP, which guarantees only solution non-negativity. A number of physical quantities (concentration, temperature, etc.) are non-negative by their nature and their approximations should be non-negative as well. We present a nonlinear finite volume (FV) method on conformal polyhedral meshes that satisfies the monotonicity condition for a wide range of problem coefficients. We admit a jumping diffusion coefficient represented by full anisotropic tensors, a jumping advection coefficient, which may be produced by the Darcy equation in multimaterial media, and both diffusion-dominated and advection-dominated regimes. The presented method is the extension of numerical schemes [9, 26] developed for the 3D diffusion equation [9] and the 2D advection–diffusion equation with continuous coefficients [26]. The major difficulty encountered in the design of a monotone numerical scheme is suppressing unwanted spurious (non-physical) oscillations in the numerical solution. These oscillations may appear in advection-dominated problems due to internal shocks and boundary layers, and in diffusion-dominated problems in highly anisotropic media due to inappropriate approximations of the diffusive flux. In the finite element (FE) context, efficient damping of spurious oscillations in advection-dominated regimes has been developed within the streamline upwind ∗Institute of Numerical Mathematics, Russian Academy of Sciences, Moscow 119333, Russia This work has been supported in part by RFBR grants 08-01-00159, 09-01-00115, Federal Program ‘Scientific and Educational Stuff of Innovative Russia’, and grant from Upstream Research Center, ExxonMobil corp. 336 K. Nikitin and Yu. Vassilevski Petrov–Galerkin (SUPG) method [5]. However, spurious oscillations around sharp layers may still appear in the SUPG solution. Spurious oscillations at layers diminishing (SOLD) methods [11] are generalizations of SUPG, which satisfy the DMP at least in some model cases. Spurious oscillations of FE solutions in diffusiondominated regimes are caused by approximation difficulties in the case of general meshes and diffusion tensors. The theoretical analysis of DMP in the FE methods [8, 15, 32] imposes severe restrictions on the coefficients and the computational mesh. An algebraic flux correction [16, 17] is the alternative approach to the design of monotone FE methods. We note, however, that many FE methods are formally not locally conservative on the cells of the original computational mesh. The finite volume (FV) methods, in contrast, guarantee the local mass conservation by their construction. The development of new FV methods for the advection– diffusion equation has been a popular topic of research, (see [3, 4, 10, 20, 28, 36]) for the steady equation and [35] for the unsteady equation and references therein). The advective fluxes can be approximated via the upwinding approach and controlled with different slope-limiting techniques [4, 7, 23] or the introduction of artificial viscosity [2, 28]. Many advanced second-order accurate linear methods for the diffusion equation fail to satisfy the monotonicity condition [1, 24, 30]. Nonlinear methods have seemed to be the feasible approach towards monotone and second-order accurate discretization [4, 11]. Nonlinear methods have been developed for the Poisson equation [6] and for the general diffusion equation [9, 14, 18, 21, 22, 24, 27, 29, 35, 37]. Our approximation of the advective flux is the 3D extension of the 2D nonlinear method proposed in [26]. The method follows the idea of the MUSCL method [34] and uses a piecewise linear discontinuous reconstruction of the FV solution on polyhedral cells, whose slope is limited via a three-by-three matrix with nonlinear entries. More precisely, we minimize the deviation of the reconstructed linear function from the given values at selected points subject to some monotonicity constraints, which form a convex set in the space of the function gradient components. The constraints are related to those considered in [12], but differ in the set of selected points and in the norm of deviation. For the discretization of the diffusive flux we adopt the nonlinear two-point flux approximation on polyhedral meshes proposed in [9]. The original idea was proposed by Le Potier in [21] for the explicit scheme for the unsteady diffusion equation on triangular meshes. Further developments of the method [14, 24, 35, 37] extend it to a wider class of meshes and equations, but inherit the interpolation from primary unknowns defined at the mesh cells to secondary unknowns at the mesh vertices. The use of interpolation affects the accuracy of the numerical scheme, as well as the properties of nonlinear solvers. An interpolation-free nonlinear FV method on 2D meshes with polygonal cells was developed in [25]. It was extended to polyhedral meshes in [9] using physical interpolation for secondary unknowns at boundary faces and faces where the diffusion tensor jumps. The proposed FV method is exact for linear solutions. Therefore, for problems with smooth solutions, one can expect the second-order asymptotic convergence A monotone nonlinear finite volume method 337 rate, which is confirmed in our numerical experiments. The monotone properties of the discrete solution are illustrated by the numerical experiments as well. The twopoint stencil for flux approximation results in sparse matrices even on polyhedral meshes. For cubic meshes and a diagonal diffusion tensor these matrices reduce to the conventional 7-point discretization. Although the method is not interpolationfree, most of the interpolation operations are based on physical principles and thus do not affect the numerical properties of the method. The major computational overhead in nonlinear FV methods is related to two nested iterations in the solution of a nonlinear algebraic problem. The outer iteration is the Picard method, which guarantees the solution non-negativity on each iteration. The set of admissible gradients in linear reconstruction used in the discretization of the advective fluxes is chosen to guarantee the stability of Picard iterations. The inner iteration is the Krylov subspace method for solving linearized problems. The paper outline is as follows. In Section 1, we state the steady advection– diffusion problem. In Section 2, we describe the construction of discrete fluxes which form the basis of our method. In Section 3, we discuss the properties of the resulting algebraic system and present our algorithm for the generation and solution of that system. In Section 4, we present the numerical properties of the scheme using tetrahedral, hexahedral, and triangular prismatic meshes. 1. Steady-state advection–diffusion equation Let Ω be a three-dimensional polyhedral domain with the boundary Γ = ΓN ∪ ΓD where ΓD ∩ ΓN = ∅ and ΓD has a non-zero measure. We consider a model advection–diffusion problem for an unknown concentration c [see 19, 31]: div (vc−K∇c) = g in Ω c = gD on ΓD (1.1) −(K∇c) ·n = gN on ΓN where K(x) is a symmetric positive definite possibly anisotropic piecewise continuous diffusion tensor, v(x) is a piecewise continuous velocity field, g ∈ L(Ω) is a source term, n is the exterior normal vector, and gD, gN are given boundary data. It is known [19] that under the above assumptions and appropriate restrictions on gD, gN , equation (1.1) has a unique weak solution c ∈W 1 0 (Ω). We denote by Γout the outflow part of Γ where v · n > 0, and define Γin = Γ \Γout. We assume that ΓN ⊂ Γout. The sufficient conditions for the non-negativity of the solution c(x) are g(x) > 0, gD > 0, and gN 6 0. We assume that these conditions hold. From a physical viewpoint, the requirements g(x) > 0 and gN 6 0 mean that no mass can be taken out of the system. The Dirichlet boundary condition on Γout and the discontinuity in boundary data on Γin may result in parabolic boundary layers. Exponential boundary layers may appear at the part of Γout where v ·n > 0. An ideal discretization scheme must add 338 K. Nikitin and Yu. Vassilevski a numerical diffusion, which is small enough to avoid excessive smearing of the boundary layers, but sufficient to damp non-physical oscillations. 2. Monotone nonlinear FV scheme on polyhedral meshes In this section, we derive a FV scheme with a nonlinear two-point flux approximation. Let q = −K∇c+ cv denote the total flux, which satisfies the mass balance equation: div q = g inΩ. (2.1) Let T be a conformal polyhedral mesh composed of NT shape-regular cells with planar faces. We assume that each cell is a star-shaped 3D domain with respect to its barycenter, and each face is a star-shaped 2D domain with respect to the face’s barycenter. We assume that T is face-connected, i.e. it cannot be split into two meshes having no common faces. We also assume that the tensor function K(x) and the velocity field v(x) vary slightly inside each cell and div v ∈ L(Ω), div v > 0 for almost every x ∈ Ω; however K and v may jump across the mesh faces, as well as may change the direction (the principal directions for K), although the normal component of vmust be continuous on any mesh face. We denote by nT the exterior unit normal vector to ∂T and by n f the normal vector to face f fixed once and for all. On a boundary face, the vector n f is exterior. We assume that |n f | = | f | where | f | denotes the area of face f . Let NB be the number of boundary faces. By FI , FB we denote disjoint sets of interior and boundary faces. The subset FJ of FI includes the faces where K(x) or v(x) jump. The set FB is further split into subsets F D B and F N B where the Dirichlet and Neumann boundary conditions, respectively, are imposed. Alternatively, the set FB is split into subsets F out B and F in B of faces belonging to Γout and Γin, respectively. Finally, FT and ET denote the sets of the faces and edges of the polyhedron T respectively, whereas E f denotes the set of the edges of the face f . Integrating equation (2.1) over a polyhedron T and using Green’s formula we get: ∑ f∈∂T χT, f q f ·n f = ∫ T f dx, q f = 1 | f | ∫